%clear

%%

theta = [0.5 : 1. : 180.]';
phi   = [0.5 : 1. : 360.]';



%  theta = [0.25 : 0.5 : 180.]';
%  phi   = [0.25 : 0.5 : 360.]';

% theta = [104.75 : 0.25 : 129.625]';
% phi   = [0.125 : 0.25 : 359.875]';

%File_1 = 'C:\Manoj\projects\stormi\test_one_day\1x1_shell_3000Ohm_q10_17.res';

%File_1 = 'C:\Manoj\projects\stormi\test\1x1_shell_3000Ohm_q10_17.res';

%File_1 = 'C:\Manoj\projects\stormi\test\1x1_shell_10KOhm_q10_17.res';
%File_1 = 'C:\Manoj\projects\stormi\test\1x1_shell_10KOhm_q11_17.res';
%File_1 = 'C:\Manoj\projects\stormi\test\run10Kdeg1_q10_23.res'
%File_1 = 'C:\Manoj\projects\stormi\test\run3K05deg_q10_01.res';
%File_1 = 'C:\Manoj\projects\stormi\test\run3K05deg_q10_17.res'
%File_1 = 'C:\Manoj\projects\stormi\test\run3K025deg_q10_17.res';
%File_2 = 'apr00_run_2_new_PC.res';
%File_1 = 'D:\Manoj\stormi\run_100ohm-m\1x1_shell_100Ohm\1x1_shell_100Ohm_q32_23.res';
%File_1 = '200004_run_1.res';
%File_2 = '200004_run_2.res';
%File_1 = 'D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\tsunami25_g77_lx.res';
%File_1 = 'C:\Manoj\temp\hemant_x3dg_files\Ocean_apr2001_reg.res';
%File_1 = 'D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\resfiles\tsu1.res';
File_1 = '/Users/manojnair/projects/aarchive_polar_circulation/Polar1.res';


fid = fopen(File_1, 'r', 'n');

LBlock = fread(fid, 1, 'long');
NO = fread(fid, 1, 'long');
Period = fread(fid, 1, 'double'); 
np = fread(fid, 1, 'long');
nt = fread(fid, 1, 'long');
theta_first = fread(fid, 1,'double'); %ONLY LAT BAND CASE
dtheta = fread(fid, nt,'double');%ONLY LAT BAND CASE
fseek(fid, LBlock, 'bof'); % LBlock if file is written on sun

LBlock
np
nt

Hpr = fread(fid, [nt np], 'float');
Hpi = fread(fid, [nt np], 'float');
Htr = fread(fid, [nt np], 'float');
Hti = fread(fid, [nt np], 'float');
Hrr = fread(fid, [nt np], 'float');
Hri = fread(fid, [nt np], 'float');



% Epr = fread(fid, [nt np], 'float');
% Epi = fread(fid, [nt np], 'float');
% Etr = fread(fid, [nt np], 'float');
% Eti = fread(fid, [nt np], 'float');
Hx1  = -(Htr+j*Hti)*100*4.*pi;
Hy1  =  (Hpr+j*Hpi)*100*4.*pi;
Hr1  =  (Hrr+j*Hri)*100*4.*pi;

% Ex1  = -(Etr+j*Eti)*1.E+06;
% Ey1  =  (Epr+j*Epi)*1.E+06;

fclose(fid);

% fid = fopen(File_2, 'r', 'n');
%%

map = [    0.2510         0    0.2510
    0.2331         0    0.3045
    0.2151         0    0.3580

    0.1972         0    0.4115

    0.1793         0    0.4650

    0.1613         0    0.5185

    0.1434         0    0.5720

    0.1255         0    0.6255

    0.1076         0    0.6790

    0.0896         0    0.7325

    0.0717         0    0.7860

    0.0538         0    0.8395

    0.0359         0    0.8930

    0.0179         0    0.9465

         0         0    1.0000

    0.0091    0.0909    1.0000

    0.0182    0.1818    1.0000

    0.0273    0.2727    1.0000

    0.0364    0.3636    1.0000

    0.0455    0.4545    1.0000

    0.0545    0.5455    1.0000

    0.0636    0.6364    1.0000

    0.0727    0.7273    1.0000

    0.0818    0.8182    1.0000

    0.0909    0.9091    1.0000

    0.1000    1.0000    1.0000

    0.2500    1.0000    1.0000

    0.4000    1.0000    1.0000

    0.5500    1.0000    1.0000

    0.7000    1.0000    1.0000

    0.8500    1.0000    1.0000

    1.0000    1.0000    1.0000

    1.0000    1.0000    0.8333

    1.0000    1.0000    0.6667

    1.0000    1.0000    0.5000

    1.0000    1.0000    0.3333

    1.0000    1.0000    0.1667

    1.0000    1.0000         0

    1.0000    0.9333         0

    1.0000    0.8667         0

    1.0000    0.8000         0

    1.0000    0.7333         0

    1.0000    0.6667         0

    1.0000    0.6000         0

    1.0000    0.5333         0

    1.0000    0.4667         0

    1.0000    0.4000         0

    1.0000    0.3333         0

    1.0000    0.2667         0

    1.0000    0.2000         0

    1.0000    0.1333         0

    1.0000    0.0667         0

    1.0000         0         0

    1.0000         0    0.0909

    1.0000         0    0.1818

    1.0000         0    0.2727

    1.0000         0    0.3636

    1.0000         0    0.4545

    1.0000         0    0.5455

    1.0000         0    0.6364

    1.0000         0    0.7273

    1.0000         0    0.8182

    1.0000         0    0.9091

    1.0000         0    1.0000];

colormap(map);



set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', 'centimeters', 'PaperPosition', [0.2 0.2 20.5 29]); 



  load coast;



subplot(2,1,1);

%imagesc(phi, theta, real(Hr1), [-4 4] );

imagesc(phi, theta, real(Hr1));

hold on; plot([long; 360+long], 90-[lat; lat], '-k');

colorbar('h');

set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);

%title 'Ex; mV/km; run = 1; apr00';

title 'Hr; nT; run = 1; test1res';



subplot(2,1,2);

%imagesc(phi, theta, real(Hr2), [-4 4] );

imagesc(phi, theta, real(Hy2), [-4000 4000] );

hold on; plot([long; 360+long], 90-[lat; lat], '-k'); hold off

colorbar('h');

set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);

%title 'Ex; mV/km; run = 2; apr00';

title 'Hr; nT; run = 2; apr00';



%print -depsc -tiff 'Re_Hr_r1_r2_apr00.eps';



return


% 
% fid = fopen('Ex_Ey_3000.txt', 'w');
% 
% fprintf(fid, '  M2 tide; rho_lithosphere = 3000; x - north, y - east  \n');
% 
% fprintf(fid, '\n');
% 
% fprintf(fid, ' theta, phi, ReEx, ImEx, ReEy, ImEy \n');
% 
% for iy = 1 : nt
% 
%     for ix = 1 : np
% 
%         fprintf(fid, '%14.5e    ',  theta(iy));
% 
%         fprintf(fid, '%14.5e    ',    phi(ix));
% 
%         fprintf(fid, '%14.5e    ',  real(Ex1(iy,ix)));
% 
%         fprintf(fid, '%14.5e    ',  imag(Ex1(iy,ix)));
% 
%         fprintf(fid, '%14.5e    ',  real(Ey1(iy,ix)));
% 
%         fprintf(fid, '%14.5e    ',  imag(Ey1(iy,ix)));
% 
%         fprintf(fid, '\n');
% 
%     end
% 
% end
% 
% fclose(fid);

% 

% fid = fopen('Ex_Ey_300.txt', 'w');

% fprintf(fid, '  M2 tide; rho_lithosphere = 300; x - north, y - east  \n');

% fprintf(fid, '\n');

% fprintf(fid, ' theta, phi, ReEx, ImEx, ReEy, ImEy \n');

% for iy = 1 : nt

%     for ix = 1 : np

%         fprintf(fid, '%14.5e    ',  theta(iy));

%         fprintf(fid, '%14.5e    ',    phi(ix));

%         fprintf(fid, '%14.5e    ',  real(Ex2(iy,ix)));

%         fprintf(fid, '%14.5e    ',  imag(Ex2(iy,ix)));

%         fprintf(fid, '%14.5e    ',  real(Ey2(iy,ix)));

%         fprintf(fid, '%14.5e    ',  imag(Ey2(iy,ix)));

%         fprintf(fid, '\n');

%     end

% end

% fclose(fid);

